The effect of phenotypic selection on stochastic gene expression 
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Genetically identical cells in the same population can take on phenotypically variable states, 
leading to differentiated responses to external signals, such as nutrients and drug-induced stress. 
Many models and experiments have focused on a description based on discrete phenotypic states. 
Here we consider the effects of selection acting on a single trait, which we explicitly link to the 
variable number of proteins expressed by a gene. Considering different regulatory models for the 
gene under selection, we calculate the steady-state distribution of expression levels and show how the 
population adapts its expression to enhance its fitness. We quantitatively relate the overall fitness 
of the population to the heritability of expression levels, and their diversity within the population. 
We show how selection can increase or decrease the variability in the population, alter the stability 
of bimodal states, and impact the switching rates between metastable attractors. 



I. INTRODUCTION 



Within one population, individual organisms often dis- 
play a large amount of observed diversity. In naturally 
occuring populations, some of the diversity is explained 
by genetic differences between the organisms. However, 
even in genetically identical populations, such as bac- 
teria or yeast grown in the laboratory [1], we observe 
phenotypic diversity, such as the variable protein levels 
in particular cells of the same population cultured in the 
same environment. This phenotypic diversity is linked 
to intrinsic molecular noise in gene expression stemming 
from relatively small copy numbers of transcription fac- 
tors and the probabilistic nature of chemical reactions. 
While molecular noise is unavoidable, imposing physical 
limits to the precision of biochemical regulatory systems, 
it may also have a functional role (2^. In particular, it 
leads to a natural diversification of a genetically identical 
and otherwise homogenous population. Such cell-to-cell 
variability can be useful for surviving in an unexpect- 
edly changing environment or large random fluctuations 
in external signals. Such arguments have been brought 
forward to explain the larger variable duration of com- 
petence in the native circuit of B. subtilis than in the 
less noisy "synex" system [31 H] . Another classical exam- 
ple is antibiotic resistance, when a fraction of bacterial 
cells become dormant by entering an antibiotic-resistant 
state without external signals, allowing the population to 
explore two different strategies [5H7j. 

Phenotypic selection under fluctuating environments 
has recently been studied theoretically |MT^. These 
studies have formalized the observation that it is benefi- 
cial for populations to "hedge their bets" against possible 
environmental stresses by keeping small, specialized sub- 
populations able to survive in various stress conditions, 
at the cost of a lower fitness in normal conditions. To 
achieve this, cells switch stochastically between different 
phenotypic states, with rates adapted to the statistics 
of environmental changes. In this description however. 



the phenotypic space is usual reduced to a discrete set of 
states, and does not account for the molecular basis of 
noise. 

Phenotypic differences can be directly linked to the 
noisy molecular nature of regulatory circuits. For exam- 
ple, in the competent system, small comK copy numbers 
are responsible for the observed noisy duration of the 
competent state [1]. The large variability of gene expres- 
sion is genetically encoded in the design of the circuit, 
for example in networks exhibiting bimodal expression 
[T]. Phenotypic variability may also take the form of 
"epigenetic" modifications, in particular on chromatin, 
which play an important role in eukaryotic cells. Unlike 
genetic variations, these different sources of phenotypic 
variability are not transmitted to the daughter cells in 
a hardwired manner. They allow populations to recover 
from environmental stress on much faster timescales than 
traditional genetic changes. As such, they allow cells to 
try out faster and more easily reversible strategies than 
genetic evolution. 

The variability of protein copy numbers in monoclonal 
populations has been extensively studied both theoreti- 
cally and experimentally |13H17j . In this paper we want 
to examine how selection acting on a population of genet- 
ically indentical, but phenotypically variable cells, shapes 
the observed variability and the stability of the pheno- 
typic states in this population. As is often done in ex- 
periments, we associate our phenotypic state with the 
protein copy numbers of a given type of protein. We ex- 
plicitly model the stochastic dynamics of this protein, the 
expression of which is under the control of a simple gene 
regulatory network. Selection acts on the population of 
cells as a function of the numbers of copies of this protein. 
This model allows us to study the effects of phenotypic 
selection on observable traits in monoclonal populations. 

We first study the effects of various types of selec- 
tion pressures on the observed distribution in a simple 
model of constitutive (unregulated) gene expression. We 
then consider a self-activating gene, which can result in a 
bistable dynamical system, and we study the effect of se- 
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FIG. 1. Model of phenotypic selection by a single gene. In 
each cell, the number of protein n undergoes a birth-death 
process describing the synthesis (with rate b) and degradation 
(with rate d) of proteins. With rate s{n), cells divide. In 
order to keep the population size constant, the new offspring 
displaces another cell picked at random, creating an global 
and uniform death rate {s{n)}. 



lection on the steady-state occupancy of the two states, 
as well as the switching rates between them. We also 
look at the effects of selection on a gene whose expres- 
sion state changes on slow timescales compared to the 
timescale for protein change, which results in a bimodal 
distribution of protein copy numbers. 



II. MODEL OF PHENOTYPIC SELECTION 

We assume that selection acts on a single trait — the 
concentration or number of copies of a given protein in 
the cell — denoted by n. The individual fitness of cells 
is defined by the n-dependent growth rate s{n). Within 
each cell, we consider the explicit dynamics of the gene 
expression network that produces the proteins governing 
the fitness of the cell. For simplicity of exposition we 
first assume that the gene producing these proteins is 
constitutively expressed. We reason directly at the level 
of proteins by assuming that the dynamics of mRNAs 
is fast. The generalization of our framework to more 
complicated modes of gene expression is straightforward, 
and we will later go beyond constitutive expression to 
model self-regulation. 

We describe the population by the mean number of 
cells pn expressing n proteins, ignoring fluctuations stem- 
ming from small numbers of cells. We will show below 
that this approximation works well as soon as the popu- 
lation is large enough. 

The change in p^ is described by a simple birth-death 
process accounting for the synthesis and degradation of 
protein molecules in each cell, and a growth rate s„ ex- 
perienced by each cell: 

dtPn = bpn-i + d{n + l)pn+i - {b + dn)pn + s(n)p^l) 
dtp = Cp (2) 



The growth rate is the net effect of cell division and cell 
death, and may be negative. Fig. [l] summarizes the pro- 
cesses governing the internal dynamics of cells as well as 
the population dynamics. More complex modes of reg- 
ulation or gene expression dynamics can be modeled by 
chosing different forms for £. 

Given the dynamics of the population in Eq. [2J the 
normalized probability of finding a cell with n protein 
copies, is given by, p„ = Pn/J2n' /^«'' ^'^d follows: 



dtp = {C - {s)l)p, 



(3) 



where (s) = X^n '^('^)p('^)- '^^^ addition of the selection 
term breaks detailed balance and introduces a nonlinear- 
ity in the master equation. A general closed-form analyt- 
ical solution cannot be found for the steady state distri- 
bution. Assuming we know the value of (s), which must 
be expressed in terms of the parameters of the problem, 
we can still write the steady state solution in the form 
of a series, because the problem is one dimensional. In 
practice we can easily find solutions numerically by iter- 
ative Euler integration. However, in certain special cases 
that we present below, we can find an analytical solution 
for the steady state distribution given (s). 

When the number of expressed proteins is large, it is 
useful to turn to a continuous description where the pro- 
tein concentration is described by a a continuous variable 
X. Expanding Eq. [3] to second order, we get for the evo- 
lution of the probability density function P{x): 



dtP{x, t)=-d^ [f{x)P{x, t)] + dl [D{x)P{x, t)] 
+ {s{x)-{s))P{x,t), 



(4) 



where f{x) = b — dx and D{x) = [h + dx)/2 are the 
effective "drift" and diffusion coefficient, respectively. 
(s) — J dxs{x)P{x,t) is the average fitness in the pop- 
ulation. f{x) and D{x) can take more general forms to 
account e.g. for self-regulation. 



III. LINEAR SELECTION 

We first consider an exactly solvable model where the 
selection pressure is linearly proportional to the number 
of protein copies in the cells, s{n) =^ sq + sn. The evolu- 
tion of the mean number of proteins is given by: 



d{n) 
dt 



b-d{n)+si{n^)-{n)^), 



(5) 



which combines the deterministic effect of birth and 
death with Fisher's relation jTS]. 

The steady state solution of Eq. |3]can readily be found 
in generating function space (see Appendix [A|for details). 
Formally, we find an infinite family of solutions for each 
possible (n) , only one of which is numerically stable (sta- 
bility is checked by evolving Eq. l3] iteratively using Eu- 



ler's integration method) . This solution is a Poisson dis- 
tribution with a rescaled mean: 



(6) 



b = 10d, s = 0.05d 



b = 10d, s = 0.3d 



At long timescales, positive selection s > acts as an 
effective anti-degradation term — it helps cells with large 
protein copy counts to survive, and eliminates cells with 
low copy numbers. As a result the mean of the Pois- 
son distribution is shifted to higher protein copy num- 
bers. Negative selection s < has the exact opposite 
effect. The impact of selection on the average fitness of 
the population is 



(s) = So + s(n) = So + 



bs 



bs^ 



(7) 



si 



d{d-s)' 



where si = Sp -|- bs/d is the mean value of s{n) when no 
selection is present. The fitness improvement scales like 
s^ > — whether selection is positive or negative, the 
population adapts to find a better place in phenotypic 
space. 

Based on this simple model, we see the general effect 
of selection that will come back in more complex sys- 
tems. If we consider the potential landscape of the reg- 
ulatory network, the system reaches a balance between 
the selection force s that is perturbing the protein con- 
centration in the cell, and the restoring force coefficient 
d due to the birth-death process. For this reason, to see 
visible and non-trivial effects of selection, the timescales 
of selection and the restoring force must be comparable. 
For very strong selection, the mean number of protein 
grows uncontrollably ( (n) -^ oo when s — >■ d) as selection 
amplifies very rare cells with abnormally large protein 
numbers. As we shall see, these effects have more visible 
consequences when regulation, and even more bimodal- 
ity, come into play. 

As a general test of our mean-field approximation, 
whereby we reduce the system to a density function p„, 
we verify our analytic result against Gillespie simulations 
of populations of cells. We explicitly consider N cells, 
in which the gene regulatory network is modeled by a 
standard time varying Monte Carlo (Gillespie) algorithm 
P^ [5D] , which appropriately models the regulation func- 
tion for the different systems we consider (constitutive 
expression, self-activation). We assume that all cells di- 
vide stochastically with rate s{n). In order to sample 
the steady-state distribution, we keep the population size 
constant, by compensating each division by the removal 
of a random cell. Fig. [2] shows the difference between the 
analytic solution and the results of the simulation for in- 
creasing population size TV, as measured by the KuUback- 
Leibler divergence (DKL). For small populations the ef- 
fect of selection is moderate, and even absent in the ex- 
treme case A^ = 1 where selection is irrelevent. As the 
population gets larger the theoretical prediction becomes 
more and more accurate. 
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FIG. 2. Validity of the description by a density function. 
Gillespie simulations of all cells in the population are com- 
pared to the analytic prediction for p„ (the fraction of cells 
with n proteins) under linear selection {s(n) = so + sn) using 
the KuUback-Leibler divergence (DKL) between probability 
distributions. Numerical results show excellent agreement for 
large population sizes. When selection is strong (s = 0.3d, 
right panel) larger population sizes are needed to reach good 
agreement than when selection is weak (s = 0.05d, left panel). 
The insets show the unselected (dashed line) and selected (full 
line) distributions of protein numbers. 



IV. SELF REGULATING GENE 

We now turn to study the effect of regulation on phe- 
notypic selection. Regulation is modeled in the simplest 
manner by assuming that the birth rate locally depends 
linearly on n with coefficient 6i = db/dn, b{n) « b + bin. 
Let us first examine the behaviour when no selection is 
present. In this case the Eq. [3] can be solved using the 
generating function technique, and the solution reads (see 
Appendix IbI) . 



= 1- 



fc/bl 



1 






(8) 



Compared to the case with no regulation, the distribution 
is no longer Poisson: the mean shifts to (n) — b/{d~bi), 
and the Fano factor is larger: ((n^) — (ri)^)/(n) = 1 + 
, ''\ . Self-regulation increases the mean and the rela- 
five variance, while self-repression decreases them both. 
Solving Eq. |3] in the presence of selection, we find again 
an infinite family of solutions. Numerical simulations 
show that the only stable solution is the one that cancels 
one of the poles of the generating function. This solution 
takes on the same functional form as Eq. [8j where d is 
replaced by a rescaled death rate defined as: 



d = 



i (d - s + 6i + v/(rf + 6i-s)'-46id) 



(9) 



which simplifies to d = d — s/(l — bi/d) in the limit of 
small selection coefficient s, and to d— s — bis/{d— s) in 
the limit of small 6i. The relative effect of selection on 



(n) can be evaluated for small s: 



(d-hy 



(10) 



and the population fitness improvement reads « 

The effect of positive regulation is to lower the effective 
restoring force to the mean value, increasing fluctuations 
in the protein copy number, as indicated by the increased 
Fano factor. These large fluctuations allow cells to ex- 
plore and find regions of larger fitness, increasing the 
mean fitness of the population. Negative regulation has 
the opposite effect. 

To gain further insight into this as well as other, more 
general models of regulation, we consider the continu- 
ous limit of the model, for which we can find an ana- 
lytic solution in the small noise approximation. Fluc- 
tuations around the steady-state value are assumed to 
be small. The mean steady-state concentration Xq is 
defined by /{xq) = 0. In the vicinity of xq, we can 
expand at leading order in the limit of small fiuctua- 
tions: D{x) ^ D{xo) = D, /{xq) — —k{x — xq) and 
s{x) w So + s{x — Xq). Then Eq. Hsimplifies to: 

dtP = kd,^[{x ~ xo)P] + DdlP + s{x - {x))P. (11) 



The steady-state solution to dtP{x, t) 
(see Appendix [A|: 



is given by 



P{x) 



y/2TrD/k 



exp 



k 
2D 



Xq 



Ds 



(12) 



As with the discrete birth-death process, the effect of 
selection is to change the mean concentration. This shift 
is proportional to the selection coefficient s, and the noise 
D. The mean population growth rate is also affected by 
this shift in a quadratic manner: 



(s) ^so+Ds^/k^ 



(13) 



The parameter k may physically be interpreted as the 
stiffness of a spring. The larger the stiffness, the less 
cells are allowed to explore regions of potentially higher 
fitness, and the smaller the advantage confered by selec- 
tion to the population. Within our small noise expan- 
sion, we have b{xQ) — dxQ, k Ki d— b'{xo) and D « b{xo)- 
The parameter fe'(xo) quantifies regulation and is equiv- 
alent to bi in the discrete birth death process. Activa- 
tion (b'{xQ) > 0) favors fiuctuations away from the mean 
steady-state value, while repression b'{xo) < suppresses 
them. The critical point b'(xQ) — d, where everything di- 
verges, marks the transition towards a bistable system, 
which we discuss in Sec. IVIl 

The scaling of the population fitness improvement 
(Eq. 13) with D and k can be interpreted as follows. 



D /k is the variance of protein number fluctuations, and 
thus quantifies the extend to which cells are allowed to 
explore better regions of the phenotypic space. k~^ is 



the relaxation time of gene expression, and quantifies 
how long cells keep the memory of their internal state, 
and how reliably they can transmit it to their offspring 
across generations, i.e. their memory of heritability. Not 
only is it important to hedge one's bets to adapt quickly 
to environmental changes, but cells must also transmit 
these fluctuations to offspring for the population to ben- 
eflt from them in the long run. The fitness improvement 
due to selection is thus the product of these two features, 
variability and heritability. 



V. THRESHOLD AND CLIFF SELECTION 

Bacterial cells grown in the presence of an antibiotic 
can develop resistance to the drug without changing the 
genome [H |5l [2TJ [22] , by expressing an antibiotic resis- 
tance gene above a given threshold. To describe this 
situation we assume that cells with at least ric of pro- 
tein copies reproduce with rate sq in the presence of the 
drug, whereas cells with n < ric grow with rate si. The 
selective pressure now takes the form of a step function, 
s{n) = Si + {sq — si)8(n — ric). We call this scenario 
threshold selection. Note that because of normalisation, 
the distribution of protein levels in the population does 
not depend on the absolute scale of s{n), and the only 
relevant parameter is As = sq — si. In the extreme case 
where cells under the threshold die, we have As = -l-oo. 
We call this scenario cliff selection. 

By inspecting the effects of threshold selection on a 
constitutively expressed gene with a mean expression 
of (n) — 20 protein copies in the absence of selection 
(Fig. [3K) , we see that even moderate selection pressures 
acting within the variance of the mean of the distribution 
result in a steep increase in the mean number of proteins. 
The cells that express small numbers of proteins now have 
a fitness disadvantage and hence the distribution shifts 
to have a higher fraction of cells expressing more pro- 
teins. For even larger selection pressures, the cells with 
n > Tic are favoured, but the mean production rate in the 
cells remains the same. Therefore a balance is reached 
between the restoring force due to protein degradation, 
which brings protein copy counts in cells down below the 
threshold to n < ric hindering their reproduction, and 
the proliferation of cells with n > ric. The mean number 
of proteins in a population thus reaches a plateau for the 
cliff model. As — >■ -l-oo. 

When As is large, as the mean number of proteins 
increases with selection pressure, the variance decreases 
and the distribution becomes subpoissonian as shown by 
the decrease in the Fano factor (see Fig. [sfe). The vari- 
ability of the population is thus reduced. Cells that sur- 
vive the selection pressure have more offspring and effec- 
tively transmit information about their expression state 
to the next generation, while cells that produce less than 
the threshold are less likely to have offspring. However 
for a relatively large critical value of ric (dash-dotted line 
in Fig. [3J3), the Fano factor increases for small selection 
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FIG. 3. The effects of threshold selection pressures on the mean (A) and Fano factor, CT^/(n), (B) of a population of cells 
expressing a constitutively expressed gene. The threshold regulation function is s = AsO(n — ric). Different values of the 
position of the threshold, ric, are compared to a system with no selection for a gene with b — 20, d = 1. Panel C shows 
examples of the distributions compared to the Poisson distribution with (n) = b/d that describes the system with no selection 
for different selection pressures for ric = 30. The variance of the distribution increases for small selection pressures and then 
the mean shifts to higher values, resulting in the initial increase and then decrease of the Fano factor. 



pressures. In this case a small fraction of cells are to 
the right of the threshold and bear a selective advantage. 
When this advantage is still small, this causes the tail 
of the distribution to get slightly fatter, thus widening 
the distribution, but without significantly affecting the 
mean. For larger selection pressures, the advantage of 
expressing more proteins becomes significant and we ob- 
serve a cusp in the probability distribution at n — ric 
(Fig. HP). 

In the limit case of cUff selection As — > +oo, where 
the effect of drugs is most detrimental, one can solve 
formally for steady state via the generating function (see 

AppendixO). As before we find a family of solutions pn , 
parametrized by a single number /3' = ricPn^ ■ /?' must be 
smaller than some critical /?, defined such that for f3' > /3, 
pii becomes negative, making the solution unphysical. 
Numerical stability analysis shows that the only stable 
solution is in fact found at this critical value /3. The 
average growth rate in the population can be calculated 
from Eq. [3] and its value is sq — (id: the proliferation 
of cells So, minus the flux of cells falling off the cliff, 
fid. Therefore /3 is the rate of cell death in units of the 
degradation rate. It is shown as a function of b/d for 
several values of Uc in Fig. [4J 



The value /3 = so/d marks the transition between 
the two phases of population: extinction and prolifer- 
ation. When P > So/d (extinction), the lifespan of the 
population under stress is given by log{No)/{Pd — sq), 
where Nq is the initial population size. Note that bi- 
ologically, in the absence of regulation, the death rate 
should be larger than the division rate because of dilu- 
tion, so/d < I. (3 = 1 therefore represents a best case 
scenario where degradation is kept to a minimum, and 
survival is maximum. The average protein level is given 
by (n) = [b/d—(3{nc—l)]/[l~(3]. Therefore the transition 
at /3 = 1 is obtained at b/d = ric— 1- 

We have seen that the effect of regulation was to rescale 
d to d — db/dn, making it possible to have so/d > I. In 
that case, the transition between extinction and prolif- 
eration is reached at higher /3, and therefore at smaller 
mean expression levels b/d (see Fig. l4J bottom left). In 
other words, positive regulation and the concomittent in- 
creased variability allow the population to better survive 
an acute stress. 

The continuous couterpart of the cliff model can also 
be solved, with f{x) = —kx, D{x) = D, and s = sq for 
X > Xc, and — oo otherwise. As in the discrete case, there 
exists a /3c above which P{x) becomes unphysical. Nu- 
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FIG. 4. Effect of cliff selection on the population. Bot- 
tom, from left to right: cell death rate, mean protein level, 
and Fano factor (variance over the mean) as a function of 
the mean unselected expression level h/d, for various values 
of Tic- Selection increases the mean protein level above ric- 
For high thresholds, cells cluster around Uc, resulting in low 
variances as shown by the Fano factor. This effect becomes 
smaller when the unselected expression is large compared to 
the threshold, b/d 2> Wc. In the leftmost plot, for a given 
proliferation rate sq, /3 = so/d marks the transition between 
extinction and proliferation. The lines can thus be inter- 
preted as separatrices between these two phases in the space 
(b/d, so/d) for various values of ric. For /? > 1 (dashed line), 
regulation is needed to achieve values of so/d that ensure 
survival. Top: example distributions of protein numbers p„ 
for the values of unselected expression levels b/d marked by 
dotted line in the bottom plots. The dashed line shows the 
location of the threshold. 



merical experiments show that this /3c is the only stable 
solution. The average growth rate of the population is 
So — PD. P = Sq/D gives the boundary in phase space 
between extinction and proliferation. The average con- 



centration is given by (a;) = (1 — k/f3D) 



indicating 



that /? = k/D when Xc = 0. At that particular point the 
solution is simply P{x) ~ kx/D exp{—kx^ /2D). 



VI. MULTISTABILITY 

In Sec. |IV| we have discussed the importance of the 
heritability of the expressed number of proteins for the 
population to benefit from selection. One of the mecha- 
nisms that has been proposed [3 [53] to stabilize pheno- 
typic states of cells with higher fitness is self-activation 
of genes. In a large parameter regime self-activating gene 
circuits are bistable. There are two deterministic steady 
state expression states: one with a high number of pro- 
tein copies and one with a low number of protein copies. 
Self-activation stabilizes these two states and leads to two 
stable subpopulations, allowing the population of cells to 
respond to different pressures. This simple scenario has 



been studied extensively in the literature [HI M HIl 121] • 
Here we consider the effects of selection on the diver- 
sity of the responses of the population, and the stability 
of each of these states, within a concrete model of gene 
expression that displays bistability through a steep self- 
regulating function, b{n) = {boK^ + biv?){K'^ + n^). 

Within the models of selection we have discussed so far, 
the high protein number state is favoured by selection. In 
Fig.jslwe show the effects of threshold (Fig.jSlA) and lin- 
ear (Fig. ^ B) selection on the mean number of proteins 
in a population of cells with bistable genes. These genes 
have a close-to-equal probability of expressing proteins 
in high and low numbers in the absence of selection. As 
discussed earlier, cells that express low protein copy num- 
bers are less likely to reproduce when selection is present. 
However bistability greatly amplifies this difference. For 
positive selection, the low protein copy expression state 
is virtually eliminated and the system looses its bistable 
nature, and hence its diversity. Analogously, negative se- 
lection pressures eliminate the high protein copy number 
expression state. This is illustrated by the probability 
distributions of protein expression shown in Fig. [Sj 

Unlike in the case of the const it utively expressed (un- 
regulated) gene, where the effects of linear selection pres- 
sure were quite smooth ((n) — b/{d — s)), the bistable ex- 
pression results in a population that is very susceptible 
to selection and a steep transition in the mean number of 
expressed proteins. Selection effectively acts on the ex- 
pression states (low/high) and does not discriminate cells 
that differ by a few numbers of proteins, resulting in a 
threshold response for linear selection. For this reason, 
the behaviour is not much affected by the precise form 
of the selection function. For large, linear selection pres- 
sures, the distribution becomes unimodal and then we 
recover the same behaviour as in the unregulated gene 
discussed in section HIH — an increase in the mean number 
of proteins as (n) = b/{d— s). In the case of threshold 
selection, for large selection pressures (positive or nega- 
tive) the system also behaves effectively like an unregu- 
lated gene, and the mean number of proteins reaches a 
plateau. 

Large threshold values Uc stabilize the expression in 
the low state, as illustrated by the dashed lines of Fig.jSK. 
Similarly to the case of no regulation discussed in Sec.|V[ 
the fraction of cells above ric is low for moderate As, 
resulting in a less fit but more diverse population than 
for lower values of Uc- 

As the mean number of proteins expressed by the genes 
increases, the response of the system to selection becomes 
steeper. The mean number of proteins expressed in the 
population in Fig. [sb is roughly double that expressed in 
the population in Fig. [5j\. The noise in the latter system 
is higher, resulting in more frequent switching between 
the low and high expression states, which can be seen 
by comparing the height of the barrier between the two 
states in the probability distributions in the absence of 
selection. As a result, low noise further amplifies the 
effects of selection by freezing the expression state in the 



lifetime of a cell, thus increasing the heritability of its 
state. 

In the limit of small noise, the system can thus be 
reduced to two states: low or high expression. If we 
know the transition rates /c+ and k- from low to high 
and from high to low, as well as the average selection 

coefficient S_ = J2n<no ^i'^)Pn and S+ = T,n>no ^i^)Pn 

in the low and high states (where no is the midpoint 
between the two states), we can write coupled equations 
for the number of cells in each of the two states, p+ = 

J2n>noPn and p_ = J2n<noPn- 



dp+ 

dt 

dp- 

~dr 



= k+p^ - k^p+ + S+P+, (14) 

= k^p+ - k+p^ + s^p^. (15) 



These equations are commonly used to describe growing 
populations with two states [8] . They have been proposed 
in the context of bacterial persistence |^ to model the 
switching between normal and persister cells in E. coli, or 
betwen the low and high expression states of a antibiotic 
resistance gene in S. cerevisiae [H], and has also been 
used in the context of the galactose utilization network of 
S. cerevisiae [23]. These equations can be readily solved 
at steady state, yielding the fraction of cells in the high 
state. 



P+ 



P+ 



vw 



As)2 + fc+As - fc + A.S 



P- 



2As 



(16) 



where k — k^ + fc_ and As — s+ — s^. When selection 
is negligible compared to the switching rate, As ^ fc, 
one recovers the equilibrium occupancy of a two-state 
model: p+ = fc_|_/fc. In the opposite limit, k <C As, where 
switching is rare, cells that are in the most favorable of 
the two states will proliferate and outcompete cells from 
the other state, and will do so much faster than they 
switch between the two states: p+ — (1/2)(1 + sign(As). 
This describes well the situation shown in Fig. [5] cells in 
the bistable population lose their diversity and all express 
high (low) numbers of proteins when selection is positive 
(negative) . 



VII. SWITCHING RATE BETWEEN 

METASTABLE STATES 

We have seen that selection could destabilise 
metastable states, especially when switching is very rare 
compared to the differences in growth rate. In that case, 
if we assume that the whole population is prepared in 
the state of lowest fitness, it typically takes only one cell 
to make the transition, in order for the whole popula- 
tion to follow suit and switch. Once that first cell has 
switched, it proliferates and its offspring quickly outcom- 
pete the cells that have remained in the state of lower 
fitness. When switching is even so rare that it is un- 
likely for a single cell out of a very large population to 



switch, Nk+ <^ 1, selection could have another, more 
subtle effect on the switching rate itself, by enhancing 
(or suppressing) the rare trajectories in gene expression 
space that make the transition. Cells that explore rare 
events towards the separatrix between the two states may 
be rewarded (or punished) by being allowed to reproduce 
(or made to die), therefore increasing (or decreasing) the 
future chance for a cell or its offspring to make the tran- 
sition. 

In practical terms, we would like to calculate the prob- 
ability that a single cell, or one of its offspring, escape 
the basin of attraction of a given state. This is a slightly 
different problem than the one we are faced with when 
dealing with a homogenous population of cells that is not 
under selection, because selection breaks detailed balance 
and favours some cells over others. Because of this, tra- 
ditional mean first passage methods are not applicable. 
However we can calculate these rates by solving Eq. [3] 
conditioned on cells not switching, which is implemented 
by a reflecting boundary condition at the midpoint no 
between the two states. By computing the rate of cells 
that would go through rig, we obtain a numerical esti- 
mate for the rate of first passage of a single cell. Fig. [6] 
shows the rates between the low and high states in the 
self- regulating bistable gene discussed in Fig. [5j3, as a 
function of the linear selection coefficient s. The effect of 
selection is to enhance transitions from the unfavorable 
to the favorable state by giving a selective advantage to 
cells that venture towards the transition point. 

To better understand this enhancement, we first con- 
sider a simplified version of the problem, where there are 
only three effective states: low, high, and an intermedi- 
ate state I0W2 between the low state and the transition 
point between the high and low states (see Fig. [7| . The 
transition rate from low to I0W2 is vfc, and from I0W2 
to high is -s/fc. Time is rescaled so that the transition 
rate from I0W2 to low is set to 1. The selective advan- 
tage (or disadvantage) along the reaction path is modeled 
by setting the growth rate to in the low state, and s 
in the I0W2 state. The population maintains a constant 
population size N. The transition rates are very low, so 
that Vk ^ 1 and kN ^ 1. Starting with all cells in the 
low state, we ask how long it take for at least one cell 
to transition to high. Before the transition happens, the 
system is described by the number of cells in the low and 
I0W2 states, Ni and N2, with N = N1+N2. We treat aU 
states where at least one cell made it to the high state 
as one big absorbing state. After some relaxation time, 
the rate of escape into the absorbing state is given by 
vfc(n2), where (712) is the average num ber of cells in the 
I0W2 state at quasi-equilibrium (c/. Eq. 16 with A:-|_ = Vk 
and fc_ — 1), and is given by VkN/{l 



3 with /c-f 

s). The rate of 



passage of the first cell to the high state is given by: 

kN 



(17) 



As s — >■ 1, cells in state I0W2 reproduce almost as fast as 
they switch back to low, providing an increasing chance 
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FIG. 5. The effects of selection on a population of self-activating bistable genes in the case of threshold {s{n) — si + AsQ(n — nc), 
panel A) and linear (s(n) = so + sn, panel B) selection. The mean number of proteins is shown as a function of the selection 
coefficient. Probability distributions for indicated values of the selection coefficient are plotted in the bottom of each panel. 
Panel A shows a comparison between two critical values of proteins, Uc = 25 which falls between the two expression states 
(solid line) and ric = 45 > (n). Regulation parameters were chosen to have nearly equal probability to be in the high and low 
expression states in the absence of selection, As, s = 0. In both cases fe(n) — ^"^2!'' 2" ■ For threshold regulation: 60 = 2, 
61 = 50, d — 1, K = 22.5. For linear regulation: bo = 2, 61 = 100, d — 1, K = 42. This value of K for the linear case was chosen 
to ensure slow switching between the two states. For smaller K the change in (n) as a function of the selection coefficient is 
even sharper. 



for switching to the high state. 

This first passage problem can also be studied within 
the small noise approximation. In this limit, the num- 
ber of protein copies x follows a random walk with drift 
f{x) and diffusion coefficient D{x) under a selection co- 
efficient s{x) (see Appendix Id]) . In the limit D{x) — )■ 0, 
the optimal reaction path can be calculated and satisfies: 
dx/dt = ± y7(a?)2^^4I)(2?)(s(xy^^7^' where (s) is the 
average fitness of the population in the basin of attrac- 
tion. The switching rate is given by the action of the 
optimal path, ~ exp(A). In the limit of small noise, this 
action reads: 



A = An 



dx[six)-{s)]/\f{x)\, (18) 



where Aq is the action in absence of selection. When 
going against a constant drift f{x), the enhancement of 
the rate is just proportional to the mean selective ad- 
vantage along the path. The stronger the adverse drift 
f{x), the smaller the enhancement. The rarity of switch- 
ing is typically affected by two factors: the strength of 
the adverse drift, and the distance to the transition point 
in phenotypic space. The enhancement of Eq. [18] is ex- 
pected to have a strong effect on transitions limited by 
long distances to the transition point and weak adverse 
drifts, and only a moderate effect on transitions limited 
by strong adverse drifts over short distances. This ex- 
plains the difference between the impact of selection on 
the two rates between the high and low states in Fig. [6] 
Although the two rates are comparable in the absence of 
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FIG. 6. Rates between the low and high states in the model 
of a self-regulating gene discussed in Fig.lsb, as a function of 
the linear selection coefficient. 




FIG. 7. A toy model for selection-aided switching. Cells 
transition from the low to high states via an intermediate 
state low2, in which they are allowed to reproduce with rate 

s. 



selection, the transition point uq « 19 is much closer to 
the low state {n « 2) than to the high state {n « 100), 
and therefore is less impacted by selection. 

Another interesting case is that of a constant stiffness, 
f{x) = —k(x — xo), and linear selection s{x) = sq + sx. 
For small s we have (s) ~ s(xo), and we get: 



A = Ao + 



3(a;final) - ^(Xinitial) 



(19) 



In this case the improvement in the switching rate is sim- 
ply proportional to the fitness difference between the ini- 
tial and final states. 

Taken together, these different estimates indicate that 
selective pressure has a significant (©(As)) effect on the 
rate of passage of the first cell. This is however a rather 
moderate effect compared to that on the steady-state oc- 
cupancy of the metastable states (Eq. 16 1. 



VIII. NON-ADIABATIC MODEL 

So far we have assumed that the binding and un- 
binding of any regulatory molecules occurs on very fast 



timescales compared to the timescale on which the pro- 
tein number changes. Experiments have shown that in 
the case of many systems the change of the gene ex- 
pression state p5H25] (from enhanced to basal expression 
and vice versa) can occur on timescales comparable with 
those on which the protein number changes. These types 
of models have been shown to result in a bimodal steady 
state distribution of protein numbers |29H31j , where one 
peak corresponds to protein expression when the gene is 
in the enhanced state and the other when the gene is the 
basal state. In this case, since the protein number and 
gene states change on comparable timescales, the protein 
number state can equilibrate in each of the gene states be- 
fore it changes. Although the detailed positions of these 
two peaks depend on the type of regulatory model (self- 
activation, self-repression, regulation by an external tran- 
scription factor protein), the general properties do not 
depend on the details of the binding rate. Therefore for 
simplicity of exposition we choose to present the problem 
for a gene that is regulated by an external transcription 
factor, resulting in a constant binding rate uj+. We then 
discuss the results for self-activation when the transcrip- 
tion factor protein binds as a dimer, 0;+ = hn^ /2, where 
h is the binding rate coefficient. 

Specifically, we consider the joint probability that the 
gene is in the enhanced (-I-) or basal (— ) expression 
state, and that n copies of the protein are present in the 
cell. Formally we have two density functions (p^,p^), 
and their associated normalized densities {p~,p'^) with 
X]n(Pn ~^Pn) ~ ^- ^c Can then write down the dynamics 
of this system as an extended birth death process, which 
also accounts for binding and unbinding of the activating 
protein: 

dtPn =J2^n!i^'Pn' + i^(^) ~ ^+]Pn + ^-Pn (20) 
n' 

dtPn = J2 A^n'^Pn' + [s(") - ^-]pt + ^+Pn, (21) 
n' 

where £^^'^ are the birth-death operators describing 
protein synthesis in the enhanced or basal gene expres- 
sion state, and a;_|_ and uj_ are the binding and unbinding 
rates of the transcription factor. 

When selection is linear, s(n) = sq -I- sn, an analytical 
solution to the steady state distribution can be found in 
generating function space [5D1 |3I] in terms of Whittakcr 
functions |32| . given that we know the mean number of 
protein copies in the system, similarly to the previously 
discussed systems. 

More intuition about the effects of selection can be 
gained from the fraction of cells that have genes in the 
enhanced state, 7r-|_ — X^nPni which is shown as a func- 
tion of selection in Fig. [21] for the unregulated gene and 
the self-activated gene, assuming transcription factors 
bind as dimers {lo^ ~ hn? /2) and threshold selection 
s(n) = So + As0(n — ric). An analysis of an effective two 
state system similar to the one presented in section |VI| 
(Eq. 14) can help us understand the probability for the 
gene to be expressed at an enhanced rate, 7r-|_, for the 
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selection coefficient As 



that have n protein copies and the gene is in the enhanced 
(p+) or basal state {p~) are plotted for different selection 
coefficients in the case of threshold regulation for a self- 
activating gene, assuming transcription factors bind as 
dimers {uj+ = hn'^/2). The change in the distributions 
are qualitatively similar for the constitutive gene. We 
explicitly see that strong positive selection favours the 
enhanced state. 

In summary, as in the case of abiabatic regulation dis- 
cussed in Sec. |VI[ selection destroys the one of the modes 
in bimodal systems, reducing the observed variability, 
even in the absence of regulation. This effect is expected 
to be stronger as the binding/unbinding rate is smaller, 
and is strongly amplified by positive regulation. 



IX. CONCLUSION 



FIG. 8. The effects of selection on a population of nonadi- 
abatic genes, slowly transitioning between on and off states. 
The probability for the gene to be found in the enhanced ex- 
pression state, 7r+, is shown as a function of the selection pres- 
sure for threshold selection pressures acting on constitutive 
genes, uj+ = u-/K, (dashed line) and self-activating genes 
with dimers binding, cj+ = hn^/2 (solid line). The inserts 
show examples of the probability distributions for s = —0.4 
(circle), s = (asterix) and s = 0.4 (cross) for the self- 
activating gene. The threshold is taken at Uc — 25, with 
b- = 2, b+ — 50, d = 1, a;_ — 0.5. K = 1 for the constitutive 
gene and h = u^/K, K — 7.7 for the self- activating gene. We 
note that in the adiabatic regime (cj ^ 1), 7r4. = 0.5 for all so 
for the constitutive gene, however (n) changes (see Fig. 3). 
The self-activating gene in the adiabatic regime is discussed 
in Fig. 6. 



constitutive gene. Summing Eq. 21 over the number of 
protein copies and solving for 7r+, we obtain: 



7^+ 



U++AsJ2r 



.Pn 



j_+As^„> (p+-^p„) 



(22) 



As As — > we recover the equilibrium result of the 
binding and unbinding rates, 7r+ — — r — • For large 
selection pressures compared to the binding/unbinding 



rates, 7r+ 



is given by the fraction of 



cells that have more proteins than the threshold and their 
genes are in the enhanced expression state and tends to 1 
for large As. Similarly, for negative selection pressures, 
(As < 0), 7r_ tends to 1 for large negative sq. 

This behaviour is shown in Fig. [8J We choose pa- 
rameters for which the probability of the gene to be ex- 
pressed in the enhanced and basal state is equal in the 
absence of selection. Selecting for a large number of pro- 
teins favours cells that are in the enhanced state and 
vice- versa. This effect, already visible for constitutive 
expression, is made more pronounced when feedback is 
present. Examples of distributions of the fraction of cells 



We have shown how selection acting on a simple phe- 
notypic trait such as the expression level of a gene, could 
significantly affect its mean expression level, diversity, 
and stability, to the benefit of the population of cells as 
a whole. 

The adaptation of monoclonal populations to challeng- 
ing environmental conditions, such as antibiotic stress or 
nutrition shortages, as studied experimentally in yeast 
pn IM] and E. coli W , is usually described by models 
of switching between a finite number of states. Our ap- 
proach goes beyond this coarse-grained description, and 
studies the effects of selection on the full spectrum of 
expression levels. In particular, we have characterized 
the stability and variability of expression within a single 
metastable state, within a simple model of constitutive 
expression. In this case, in the small noise approxima- 
tion, Eq. [13] quantifies how the population improves its 
overall fitness proportionally to the heritability k^^ and 
the variability D/k of fluctuations in protein copies. Her- 
itability can be enhanced by means of positive regulation, 
which decreases the relaxation rate k. 

When regulation is strong fc < 0, the system can be- 
come bistable, with two states of low and high expres- 
sion level. This is a case of very strong heritability, in 
which cells can transmit their expression state to their 
offspring over many generations. We have shown that 
selection destroys the bimodality of the distribution of 
gene expression, by favoring the state of highest fitness. 
This effect is all the more important when differences in 
growth rate between the states are large compared to the 
switching rates. 

In multistable systems, selection decreases the variabil- 
ity of a population by favoring some metastable states 
over others. However, within a single metastable state, 
a linear selection in the expression level mostly affects 
the mean and stability of expression, but not its vari- 
ance. By contrast, when selection is step-like, with a 
different growth rate below or above a given threshold of 
expression, selection may increase or decrease variability, 
depending on the strength of selection. Very stringent 
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selection tends to decrease the diversity of expression at 
the cost of fitness (Fig. [4]), while a moderate selection 
acting on the tail the distribution increase the variance 
by amplifying these tails (Fig. [s]) . 

In bistable systems, selection has another overlooked 
effect, which cannot be grasped by a simple two-state 
model: it enhances or suppresses the rate of switching be- 
tween the two states, by giving a selective (dis) advantage 
to cells going along the transition path. This selection- 
aided switching could serve as a mechanism for driving 
and stabilizing a population of cells through differentia- 
tion using a gradual selective pressure, for example dur- 
ing developement where phenotypic noise plays an im- 
portant role [33] . 

Our results show that selective pressure acting on the 
expression of a single gene may strongly affect its be- 
haviour at the population level. It would be interesting 
to test this idea experimentally, by measuring the proper- 
ties of gene expression (mean, variance, switching rates) 
in selective against non selective environments, for differ- 
ent modes and strength of regulation. For example such 
experiments could test the prediction that positive regu- 
lation enhances the effect of selection on the population 
mean. 

Our approach provides a broad framework for address- 
ing the effect of selection on observable phenotypic traits 
in genetically homogeneous populations, with straight- 
forward generalisations to arbitrary phenotypic spaces 
with multiple genes. 
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Appendix A: Solution for linear selection 

To calculate the steady state distribution for a lin- 
ear selection pressure s(n) = sn, we define the gener- 
ating function for the probability distribution, as G{z) = 
Sn ^"'Pn- In generating function space and steady state, 
assuming we know {s{n)) = s{n), Eq. plbecomes: 



{hz-b)G-(dz-d)—- =0, 
dz 



where h = b + s{n) and d — d— s. Eq. 
by direct integration to give 



Al 



G(z) 



i(z-l) 



1 



1 



(Al) 
can be solved 

(A2) 



where 5 = | and (3 = s ^~'|^"j+"<"^ Note that this ex- 
pression for the generating function self-consistently sat- 
isfies G"(l) — (n), so that (n) is not constrained by the 
condition of stationarity. We thus have a family of solu- 
tions, parametrized by (n) or equivalently by /3. An es- 
pecially simple solution is given by the condition (3 = 0, 



which yields the generating function of a Poisson distri- 
bution: 



G{z) 



j(^-i) 



(A3) 



This solution is the one we obtain by numerical integra- 
tion of Eq. [3] at steady-state. 

The Fokker-Planck equation for the continuous case 
(Eq|Tl|) is solved at steady state by going to Fourier 
space: 



P{p) = / dxe'P''P{x). 



(A4) 



Eq 11 then becomes at steady state (we set xo = with 
no loss of generality) : 

{p + is)dpP + {Dp^ + sx)P = 0, (A5) 

where x = J dxxP{x). The solution to this equation is: 

s{Ds/k^-x)/k 



Pip) 



_1 Dp^ 

e 2 k 



1 



pk 

s 



(A6) 



As in the discrete case, the only stable solution corre- 
sponds to a; = Ds/h/^. In this case, the exponent in the 
second term cancels, and we obtain the Fourier transform 
of a Gaussian distribution of mean Dsjk^ and variance 
D/k. 



Appendix B: Solution with self-regulation 

Here we give some details for the calculations of self- 
regulation (Sec. 



IV I . The birth rate is assumed to depend 



on n: b + bin. 

we get the steady-state solution in absence of selection 

G'iz) b 



Jsing the generating function technique, 

ion: 

(Bl) 



G{z) d — biz 



from which we infer: 
P{n) = ( 1 



bi\ lib 



n\\d 



nh 



j=0 



M 



(B2) 



With selection the equation for the generating function 
G{z) reads: 



G'{z) 



b + sin) — bz 



G{z) {d-biz){l- z) + sz' 
The right-hand side has two poles: 

bi+d~ s± a/(6i +d-s)2 -ibid 



z± 



2bi 



(B3) 



(B4) 



By analogy with the unregulated case, we make the hy- 
pothesis that the only stable solution is such that the 
pole at z_ disappears. This is satisfied if: 



(n) 



2b 



d- s + bi + ^/{bl +d- sy -ibid ' 



(B5) 
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Then we simply have: 

G'iz) 



reads: 



(B6) 



Giz) h{z+-zy 
so we get the same form as Eq. |B1[ after replacing d by 



d = hiz+ = -U- s + bi + ^/{d + bi- s)^ - Abid^ . 

(B7) 
We checked numerically that our assumption about the 
cancelation of the z_ pole in G' /G was correct. 



Appendix C: Solution for cliff selection 

This appendix contains details of the calculations in 
the model of cliff selection. We start with the discrete 
case, for which the evolution equation reads: 

dtPn ^b{pn^l-Pn)+d{{n+l)pn+l-npn+f3)pn, (CI) 

for n > Tic and p„ = for n < ric- (3 — ricPn^- The 
last term comes from the normalisation condition and 
compensates the loss of cells off the cliff, which happens 
with rate df5. The generating function can be calculated 
as a function of /3 at steady state: 



Gpiz) = e'^'il- zr / dy(3y 
Jo 



."c-lg-ai/(' 



il~y)-^-\ (C2) 



where a — b/d. Note that the form above automatically 
satisfies G{z) ~ p„^z"= as z — >■ 0, and G(l) = J2nPn — 1- 
Therefore /3 is unconstrained and entirely determines the 
solution. Guided by numerical simulation, we hypothe- 
size that the only stable solution corresponds to the high- 
est possible /3 that does not entail p„ < for some n. 
An analogous analytical solution exists for the threshold 
model, with an additional continuity condition between 
the two intervals (0, Uc — 1) and (jic-, -foo). 

In the continuous case, the Fokker-Planck equation 
reads: 



dtP = kd^ (xP) + DdlP + DI3P, 



(C3) 



with f3 = dxP\x=x^- The last term corresponds to the 
flux of cells crossing the threshold. The formal solution 



Ppi^) 



= xe-y^""^ 



u{yc)m{y(x)) - m{yc)u{y{x)) 
u{yc)Nm - myXc)Nu 



(C4) 



where y{x) 



kx^/2D, yc 



kxl/2D, Na = 



J dx xe y^''"'u{y{x)), with a = u or m. The func- 
tions m{x) and u{x) are defined as: m{x) — M{{k — 
/3L>)/2fc,3/2,x) and m(x) = U{{k - /3D)/2k, 3/2, x), 
where M and U are the confluent hypergeometric func- 
tions of the first and second kind, respectively. 



Appendix D: Optimal switching path 

Here we detail the calculation of the optimal reaction 
path of a stochastic process under selective pressure. We 
assume that all cells are equilibrated in one metastable 
state, and we consider the probability of rare paths out of 
this state. The probability of a path is given by the usual 
expression for the action, multiplied by a term reflecting 
the historical fitness of the cell relative to the rest of the 
population [11.^ 



[f-/(-)F 



P{{x{t)} ^ exp / "" dt I -^^\^^^ + s{x) - {s) 
exp I — / dt C 

■/tinitial 

(Dl) 
The Lagrangian C can be rewritten as: 



C 



AD{x) 



dx g{x) - f[x) 
dt 2D{x) 



(D2) 



with g{x) = ±V/(a;)2 - 4.D{x){s{x) - (s)). Note the 
second term in the integrand does not depend on the par- 
ticular path taken, and that the first term can be made 
arbitrarily small by setting dx/ dt = g{x) and by choosing 
the sign of g appropriately |341 135) . We are considering 
rare paths, which move against the drift, e.g. dx/dt > 
and f{x) < 0. Then the action of the optimal path reads: 



A 



""-' ^^ \f{x)\ + ^f{xr~ADix)isix)-{.s)) 
^Dix) 



(D3) 
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